宏基因组实战3. MEGAHIT组装拼接及quast评估
前情提要
如果您在学习本教程中存在困难，可能因为缺少背景知识，建议先阅读本系统前期文章
测试数据
刘博士帮助把测试数据建立了一个百度云同步共享文件夹，有非常多的好处，请读完下文再决定是否下载：
下载被墙的数据；很多数据存在google, amazon的部分服务器国内无法直接下载，而服务器一般科学上网不方便，下载数据困难。大家下载失败的数据请到共享目录中查找；
预下载好的软件、数据库；有很多需要下载安装、注册的软件(在线安装包除外)，其实已经在共享目录了，节约小伙伴申请、下载的时间；
数据同步更新；任何笔记或教程不可避免的有些错误、或不完善的地方，后期通过大家的测试反馈问题，我可以对教程进行改进。共享目录不建议全部下载或转存，因为文件体积非常大，而且还会更新。你转存的只是当前版本的一个备份，就不会再更新了。建议直接在链接中每次逐个下载需要的文件，也对文件有一个认识过程。
方便结果预览和跳过问题步骤；服务器Linux在不同平台和版本下，软件安装和兼容性问题还是很多的，而且用户的权限和经验也会导致某些步骤相关软件无法成功安装(有问题建议选google、再找管理员帮助；想在群里提问或联系作者务必阅读《如何优雅的提问》)。在百度云共享目录中，有每一步的运行结果，读者可以下载查看分析结果，并可基于此结果进一步分析。不要纠结于某一步无法通过，重点是了解整个流程的分析思路。
最后送上本教程使用到的所有文件同步共享文件夹链接：http://pan.baidu.com/s/1hsIjosk 密码：y0tb 。
MEGAHIT
https://2017-cicese-metagenomics.readthedocs.io/en/latest/assemble.html
进入我们的工作目录
安装程序
git clone https://github.com/voutcn/megahit.git
cd megahit
makecurl下载测序数据，或在百度云中下载，或使用在上节中K-mer trim的结果文件
cd ../data
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948.abundtrim.subset.pe.fq.gz
curl -O https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249.abundtrim.subset.pe.fq.gz开始组装
cd ..
mkdir assembly
cd assembly
ln -fs ../data/*.subset.pe.fq.gz .
../megahit/megahit --12 SRR1976948.abundtrim.subset.pe.fq.gz,SRR1977249.abundtrim.subset.pe.fq.gz -o combined测试文件为了方便演示，只取了原数据的一小部分，原作者用15min，我的服务器运行只用了4min。原始数据使用三种主流软件分析，运行所消耗时间、内存比较。
| Program | Time | Memory (GB) |
|---|---|---|
| IDBA-UD | 33h 54m | 123.84 |
| SPAdes | 16h 2m | 381.79 |
| MEGAHIT | 1h 53m | 33.41 |
查看拼接结果
less combined/final.contigs.fa评估组装结果
https://2017-cicese-metagenomics.readthedocs.io/en/latest/assembly-evaluation.html
安装评估软件quast
cd ..
git clone https://github.com/ablab/quast.git -b release_4.5
export PYTHONPATH=$(pwd)/quast/libs/运行QUEST
cd assembly
mkdir quast-evaluation
cd quast-evaluation
ln -fs ../combined/final.contigs.fa megahit.contigs.fa
../../quast/quast.py megahit.contigs.fa -o megahit-report
cat megahit-report/report.txt下载metaSPAdes结果评估并比较
curl -LO https://osf.io/h29jk/download
mv download metaspades.contigs.fa.gz
gunzip metaspades.contigs.fa.gz
../../quast/quast.py metaspades.contigs.fa -o metaspades-report
cat metaspades-report/report.txt
# look at the two reports in parallel
paste *report/report.txt结果如下：
Assembly megahit.contigs metaspades.contigs
# contigs (>= 0 bp) 7904 4112
# contigs (>= 1000 bp) 2763 1843
# contigs (>= 5000 bp) 582 583
# contigs (>= 10000 bp) 191 244
# contigs (>= 25000 bp) 18 43
# contigs (>= 50000 bp) 2 17
Total length (>= 0 bp) 13222363 12090326
Total length (>= 1000 bp) 11149439 11320830
Total length (>= 5000 bp) 5893043 7955570
Total length (>= 10000 bp) 3186708 5596677
Total length (>= 25000 bp) 663719 2500084
Total length (>= 50000 bp) 112488 1603525
# contigs 3847 2280
Largest contig 61397 261464
Total length 11895322 11615922
GC (%) 46.29 46.27
N50 4924 9303
N75 2524 3937
L50 594 266
L75 1455 754
# N's per 100 kbp 0.00 0.00结果N50和N75在metaspades结果更好，如果有计算资源，且不缺时间，推荐使用metaspades。但如果没有上T内存的服务器，项目周期又紧张，直接用metahit出结果。
猜你喜欢
热文：图表规范 DNA提取发Nature
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内五十位PI，六百多名一线科研人员加入。参与讨论，获得专业指导、问题解答，欢迎分享此文至朋友圈，并扫码加创始人好友带你入群，务必备注“姓名-单位-研究方向-职务”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决推荐生信技能树-微生物组版块(http://www.biotrainee.com/forum-88-1.html) 发贴，并转发链接入群，问题及解答方便检索，造福后人。
学习16S扩增子、宏基因组科研思路和分析技术，快关注“宏基因组”
点击阅读原文，跳转最新文章目录阅读